Global climate models violate scaling of the observed atmospheric variability 
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We test the scaling performance of seven leading global climate models by using detrended fluctua- 
tion analysis. We analyse temperature records of six representative sites around the globe simulated 
by the models, for two different scenarios: (i) with greenhouse gas forcing only and (ii) with green- 
house gas plus aerosol forcing. We find that the simulated records for both scenarios fail to reproduce 
the universal scaling behavior of the observed records, and display wide performance differences. 
The deviations from the scaling behavior are more pronounced in the first scenario, where also the 
trends are clearly overestimated. 
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Confidence in the simulation and prediction skills of 
global climate models (coupled atmosphere-ocean general 
circulation models [jl], j|] AOGCMs) is a crucial precon- 
dition for formulating climate protection policies. The 
models provide numerical solutions of the Navier-Stokes 
equations devised for simulating meso-scale to large-scale 
atmospheric and oceanic dynamics. In addition to the ex- 
plicitly resolved scales of motions, the models also con- 
tain parameterization schemes representing the so-called 
subgrid-scale processes, such as radiative transfer, turbu- 
lent mixing, boundary layer processes, cumulus convec- 
tion, precipitation, and gravity wave drag. A radiative 
transfer scheme, for example, is necessary for simulat- 
ing the role of various greenhouse gases such as CO2 and 
the effect of aerosol particles. The differences among the 
models usually lie in the selection of the numerical meth- 
ods employed, the choice of the spatial resolution ||, and 
the subgrid-scale parameters. 

Two scenarios (apart from a control run with fixed 
CO2 content) have been studied by the models, and the 
results are available from the IPCC Data Distribution 
Center [Q. In the first scenario, one considers only the 
effect of greenhouse gas forcing. The amount of green- 
house gases is taken from the observations until 1990 and 
then increased at a rate of 1% per year. In the second sce- 
nario, also the effect of aerosols (mainly sulphates) in the 
atmosphere is taken into account. Only direct sulphate 
forcing is considered; until 1990, the sulphate concentra- 
tions are taken from historical measurements, and are 
increased linearly afterwards. The effect of sulphates is 
to mitigate and partially offset the greenhouse gas warm- 
ing. Although this scenario represents an important step 
towards comprehensive climate simulation, it introduces 



new uncertainties — regarding the distributions of nat- 
ural and anthropogenic aerosols and, in particular, re- 
garding indirect effects on the radiation balance through 
cloud cover modification, etc. 

All of the models are capable, to varying extents, to 
reproduce the current mean state of the atmosphere ||. 
The models have been validated by comparing to histor- 
ical data and by intercomparison of the models |?], ||. 
The efforts have been restricted to traditional time se- 
ries analysis which generally assumes that the statistical 
properties of a signal remain the same throughout the 
entire period. This assumption of stationarity, however, 
is certainly not valid for climate records due to imposed 
effects such as global or urban warming. 

In our evaluation of the models, we apply detrended 
fluctuation analysis (DFA) || which can distinguish 
between trends and correlations and thus reveal trends as 
well as long term correlations very often masked by non- 
stationarities. Recently, Koscielny-Bunde et al. [HI [T^ ] 
applied DFA and wavelet techniques (see, e. g. [13|) to 
investigate temporal correlations in the atmospheric vari- 
ability. Considering maximum daily temperature records 
of various stations around the globe, they analyzed the 
temperature variations from their average values and 
found that the persistence, characterized by the corre- 
lation C(s) between temperature variations separated by 
s days, decays with a power law, 
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with roughly the same correlation exponent 7 = 0.7 for 
all stations considered. The range of this persistence law 
exceeds ten years, and there is no evidence for a break- 
down of the law at even larger time scales. Indications 
for the long term persistence through spectral analysis 
have also been obtained 14 15 1. Since the persistence 
scaling law appears to be universal, i, e. independent 
of the location and climatic zone of the stations, we use 
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it in the following for assessing the performance of the 
AOGCMs. 

For the test, we consider monthly averages of the 
daily maximum temperature from seven AOGCMs: 
GFDL-R15-a (Princeton), CSIRO-Mk2 (Melbourne), 
ECHAM4/OPYC3 (Hamburg), HADCM3 (Bracknell, 
UK), CGCM1 (Victoria, Canada), CCSR/NIES (Tokyo), 
NCAR PCM (Boulder, USA) (see § for details). We ex- 
tracted the data for six representative sites around the 
globe (Prague, Kasan, Seoul, Luling/Texas, Vancouver, 
and Melbourne). For each model and each scenario, we 
selected the temperature records of the 4 grid points clos- 
est to each site, and bilinearly interpolated the data to 
the location of the site. We analyze both scenarios but 
we focus more on the better established first scenario. 

We analyze for each site the variations ATj of the 
monthly temperatures Tj from the respective monthly 
mean temperature Tt that has been obtained by av- 
eraging over all years in the record. Quantitatively, 
persistence in AX^ can be characterized by the (auto) 
correlation function, C(s) = (ATiAT i+s ) = (l/(N - 
s )) J2i!=L S ATiATi+ s , where N is the total number of 
months in the record. A direct calculation of C(s) is hin- 
dered by the level of noise present in the finite tempera- 
ture series, and by possible nonstationarities in the data. 
Following Refs. || and f7), we do not calculate C(s) di- 
rectly, but instead study the fluctuations in the tempera- 
ture "profile" Y n = J27=i ATi. To this end, we divide the 
profile into (nonoverlapping) segments of length s and de- 
termine the squared fluctuations of the profile (specified 
below) in each segment. The mean square fluctuations, 
averaged over all segments of length s, are related to the 
correlation function C(s) (see below). In our test, we 
employ a hierarchy of methods that differ in the way the 
fluctuations are measured and possible nonstationarities 
are eliminated (for a detailed description of the methods 
we refer to jlO|) : 

(i) In the (standard) fluctuation analysis (FA) , we cal- 
culate the difference of the profile at both ends of 
each segment. The square of this difference repre- 
sents the square of the fluctuations in each segment. 

(ii) In the "first order detrended fluctuation analysis" 
(DFA1), we determine in each segment the best 
linear fit of the profile. The standard deviation 
of the profile from this straight line represents the 
square of the fluctuations in each segment. 

(iii) More generally, in the "nth order DFA" (DFAn), 
we determine in each segment the best nth-order 
polynomial fit of the profile. Again, the standard 
deviation of the profile from these polynomials rep- 
resents the square of the fluctuations in each seg- 
ment. 

The fluctuation function F(s) is the root mean square of 
the fluctuations in all segments. For the relevant case of 
long-term power-law correlations given by Eq. El with < 



7 < 1, the fluctuation function F(s) increases according 
to a power law , 

F(a)~a« > " = 1-|. 

For uncorrelated data (as well as for short-range corre- 
lations represented by 7 > 1 or exponentially decaying 
correlation functions), we have a = h. For long term 
correlations we have a > -z. 
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FIG. 1: Results of FA and DFA for the monthly averages of 
the daily maximum temperature of Luling (a-d) and Kasan 
(e-h). In each panel, the curves from top to bottom repre- 
sent the fluctuation function F(s) obtained from FA (A) and 
DFA1-DFA5 (o), respectively. Full symbols represent past 
data and open symbols are for the entire simulation period. 
The scale of F(s) is arbitrary. The maximum value of s is 1/4 
of the length of the considered record. The two lines shown at 
the bottom are theoretical lines with slope 0.65 (upper line) 
and 0.5 (lower line). Note that FA (which does not remove 
trends) overestimates the fluctuation exponent as can be seen 
when comparing to DFA. As seen in the figures, the differ- 
ences in the exponents obtained by DFA3 and higher orders 
of DFA are negligible, which means that in DFA3 all trends 
are removed from the data. For this reason, when focusing 
on the correlation exponent, it is sufficient to use DFA3 (as 
done in Figure ^) . 

By definition, the FA does not eliminate trends, sim- 
ilar to the Hurst method and the conventional spectral 
analysis [fill. In contrast, DFAn eliminates polynomial 
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trends of order n — 1 in the original data. Thus, from 
the comparison of the various fluctuation functions F(s) 
obtained by these methods, we can learn both about long 
term correlations and types of trends, which cannot be 
achieved by the conventional spectral analysis. For test- 
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FIG. 2: Results of DFA3 applied to the monthly averages of 
the daily maximum temperature of (a) real data and (b-h) 
simulated data (for scenario (i)) at the respective geograph- 
ical positions of Prague (Pr), Melbourne (Me), Luling (Lu), 
Seoul (Se), Kasan (Ka) and Vancouver (Va), for the seven 
AOGCMs considered. In each panel, curves from top to bot- 
tom represent the result for the six sites. For better evaluation 
of the models, F(s) is divided by s 0,5 . The scale of F(s) is 
arbitrary. The theoretical line at the bottom has a slope of 
0.15. 

ing the performance of the models we plot, in a dou- 
ble logarithmic presentation, F(s) versus s for FA and 
DFA1-DFA5 for each of the six sites and compare the 
curves with those obtained from the real data. Figure [j] 
shows representative results of the fluctuation functions 
for Kasan (Russia) and Luling (Texas) for both the real 
data and the data from six of the climate models with sce- 
nario (i) (three models for each location) . The data from 
the models are taken (a) until 1990 (full symbols), which 
is the simulation period corresponding to the observed 
data and (b) for the entire simulation period including 
future data (open symbols). 

Every panel shows, from top to bottom, the fluctua- 



tion functions F(s) obtained from FA and DFA1-5. We 
also have drawn two straight lines with slopes 0.5 (corre- 
sponding to uncorrelated data) and 0.65, corresponding 
to the correlation exponent obtained for the real data. 
Figure [l] shows that for the real data, for both Kasan 
and Luling, all the F(s) curves are parallel straight lines 
with slope close to a = 0.65, beyond 1 year. This indi- 
cates (i) the absence of trends and (ii) the existence of 
long term power law correlations consistent with earlier 
findings IJ. 

The simulated records show a quite different behavior. 
For Luling, CSIRO-Mk2 and ECHAM4/OPYC3 yield FA 
curves that are not parallel to the DFA curves, having a 
larger asymptotic slope, while the DFA curves show a 
crossover towards uncorrelated behavior (a — 0.5) after 
roughly 2 years. This indicates (i) an overestimation of 
the trends and (ii) the loss of long term correlations in the 
models. The HADCM3 model performs slightly better, 
with DFA curves approaching a slope of a = 0.62 at long 
times. However, when compared to real data, the FA 
curve bends slightly upwards at long times (for the past 
data) , revealing also an overestimation of the trend. For 
Kasan, the CGCM1 model yields uncorrelated behavior 
at long times. CCSR/NIES and NCAR PCM show long 
term persistence, with an exponent a slightly below 0.6 
in both cases. In all cases, the FA curves are straight 
lines, with slightly larger exponents than the DFA curves 
for CGCM1 and NCAR PCM. This again points to an 
overestimation of the trends by the models. 
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FIG. 3: Histogram of the fluctuation exponent (a) values 
obtained for the six sites simulated by the seven AOGCMs 
(a) scenario (i) and (b) scenario (ii) for the entire records. 
The entries in each box represent 'Model - Site'. 

For reviewing the scaling performance of the models 
we concentrate now on the results of DFA3, since DFA3, 
DFA4 and DFA5 yielded the same scaling exponents for 
all cases. We find that the fluctuation function in DFA3- 
5 remained unchanged while FA is changed dramatically 
when including also future data (see Figure |l|) . This fea- 
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ture shows some internal consistency of the models, since 
past and future differ mainly in the amount of trends cre- 
ated by greenhouse gases (as shown by the FA curves), 
and trends are well eliminated by higher order DFA [E0| . 
Therefore, we use the entire data for a more accurate 
estimation of scaling exponents. 

Figure |^ shows the results for the fluctuation function 
obtained from DFA3 for all available real and model data 
(all seven AOGCMs with scenario (i)) at the six sites 
considered, for the entire simulation period. To facilitate 
the evaluation of the models, we have divided F(s) by 
s 5 . A plateau now indicates loss of long term correla- 
tions. The straight line in each panel has a slope of 0.15, 
corresponding to the universal exponent in the original 
curves. While the real data yield, for long times, parallel 
lines with a slope of 0.15 for all sites in agreement with 
the earlier findings, the virtual data display wide per- 
formance differences and fail to reproduce the universal 
features of the benchmark time series. Direct inspection 
of the figure shows that about half of the model curves 
are very close to a plateau, yielding uncorrelated behav- 
ior above approximately 2 years. 

The actual long term exponents a for the greenhouse 
gas only scenario obtained by the seven models for the six 
sites are summarized in a histogram in Figure ^|a. The 
histogram shows a pronounced maximum at a = 0.5. For 
best performance, all models should have exponents a 
close to 0.65, corresponding to a peak of height 42 in the 
window between 0.62 and 0.68. Actually more than half 
of the exponents are close to 0.5, while only 6 exponents 
are in the proper window between 0.62 and 0.68. 

Figure |^b shows the histogram for scenario (ii), where 
in addition to the greenhouse gas forcing, also the effects 
of aerosols are taken into account. For this case, there is 
a pronounced maximum in the a window between 0.56 
and 0.62 (more than half of the exponents are in this 



window) , while only 5 exponents are in the proper range 
between 0.62 and 0.68. This shows that although the 
second scenario is also far from reproducing the scaling 
behavior of the real data, its overall performance is better 
than the performance of the first scenario. 

To summarize, we have presented evidence that 
AOGCMs fail to reproduce the universal scaling behav- 
ior observed in the real temperature records. Moreover, 
the models display wide differences in scaling for differ- 
ent sites. When comparing the two scenarios, our results 
suggest that the second scenario (CO2 plus aerosols) ex- 
hibits better performance regarding the values of the scal- 
ing exponents as well as the trends. The effect of aerosols 
not only decreases the trends but also modifies the fluc- 
tuations, to more closely resemble the real data. This 
confirms in a way independent of the evaluations made 
so far H that the incorporation of aerosols is necessary 
to approach reality. 

It is possible that the lack of long-term persistence is 
due to the fact that certain forcings like volcanic erup- 
tions or solar fluctuations have not been incorporated 
in the models. However, we cannot rule out that sys- 
tematic model deficiencies (such as the use of equivalent 
CO2 forcing to account for all other greenhouse gases or 
inaccurate spatial and temporal distributions of sulphate 
emissions) prevent the AOGCMs from correctly simulat- 
ing the natural variability of the atmosphere. 
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